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I. INTRODUCTION 

The phenomenon of particle localization in a disor- 
dered medium has attracted a considerable interest start- 
ing with the famous Anderson paper [1]. The scaling the- 
ory [2,3] predicts a transition from extended to localized 
electronic states for any spatial dimension d > 2 and 
sufficiently large disorder. It assumes that this transi- 
tion (Anderson transition, localization transition, metal- 
insulator transition (MIT)) can be described as a phase 
transition of second order with just one single relevant 
scaling field. Wegner [4] introduced the so-called A''- 
orbital model which allowed for a field-theoretical refor- 
mulation of the problem [5-7] followed by a renormaliza- 
tion group analysis in 2-1- e dimensions [7-9]. It was recog- 
nized that a usual mean-field approach is not possible in 
this field theory since the conventional one-point Green 
function - which is closely related to the mean density of 
states (DOS) - does not exhibit a critical behavior. This 
seems to break the Goldstone theorem. [10] By extend- 
ing and making rigorous this earlier work (which made 
use of the mathematically ill-defined replica trick) Efe- 
tov [11] mapped the problem of disordered one-electron 
systems on supersymmetric non-linear cr-models (in the 
following called supermatrix models) with certain super 
coset spaces G/K according to the particular behavior 
of the disordered one-electron system under time reversal 
and/or spin rotation. 

The investigation of the aforementioned supermatrix 
models on the Bethe lattice (BL) [12-14] revealed an un- 
expected behavior at the transition point. Namely, a 
jump in the inverse participation ratio (IPR) on the lo- 
calized side and an exponential decrease of the diffusion 
constant on the metallic side of the transition was found. 
These results gave rise to the hypothesis [14] that the 
scaling theory (which predicts a power-law behavior of 
various quantities) might be not applicable for the lo- 



calization transition. In order to investigate the critical 
behavior of this transition on a hypercubic lattice, a La- 
grangian was constructed in Ref. [15] so that its saddle- 
point reproduced the so-called effective-medium approx- 
imation (EMA) of Ref. [16] which leads to the same non- 
power-like critical behavior as in the case of the BL. It 
was proposed [17] that this exotic behavior is character- 
istic for the Anderson transition on a d-dimensional hy- 
percubic lattice for a (sufficiently high) spatial dimension 
d>2. 

However, Mirlin and Fyodorov [18] argued that con- 
trary to the above hypothesis the non-power-like critical 
behavior is an artifact of the special lattice structure of 
the BL and of the EMA which imitates this structure. 
One can assign an effective dimension d = oo to this 
lattice structure, which plays the role of the upper crit- 
ical dimension. (A transition behavior with critical ex- 
ponents and c» is then formally identified with a jump 
and an exponential behavior, resp.) In this picture the 
above critical behavior on the BL or within the EMA can 
be understood within the one-parameter scaling assump- 
tion; for finite dimensions a power-like scaling behavior 
is expected. 

Moreover, it was shown in a quantitative way [19] that 
the EMA is an uncontrolled approximation due to the 
neglect of loop graphs, since in the context of a strong- 
coupling expansion the highly branched graphs with a 
high number of loops yield the dominant contribution 
whereas zero-loop graphs can be neglected. It is there- 
fore perhaps fair to say that instead of trying to calculate 
corrections to the BL results it seems more promising 
to investigate the localization transition problem (within 
the supersymmetry formalism) directly in li = 3 dimen- 
sions, at least numerically. This has, however, not yet 
been done. 

The development of a supersymmetric toy model, the 
Hyperbolic Superplane (HSP) [20], allowed for some tech- 



1 



nical simplifications in the investigations of non-linear a- 
models. It was shown in [20] that this model is capable of 
describing localization in a quasi one-dimensional geom- 
etry. On the BL, an analysis of a non-linear cr-model 
taking the HSP as target space [19] reproduced the pre- 
viously reported unconventional critical behavior. This 
experience indicates that the HSP might yield a quali- 
tatively correct picture of the localization transition in 
d > 2 dimensions, too. 

Technical simplifications compared with the aforemen- 
tioned supermatrix spaces originate from the fact that 
the HSP has only one radial coordinate whereas the su- 
permatrix models show at least two radial degrees of free- 
dom. It may thus serve as a useful toy model which al- 
lows one to study many interesting features related to 
the MIT with relative ease, thereby gaining some expe- 
rience and insight into the main difiiculties before one 
turns to the more interesting supermatrix spaces. (The 
HSP has been successfully applied to several other prob- 
lems as well, e.g. Migdal-Kadanoff renormalization [21] 
or quantum chaos in conjunction with a superanalog of 
the Selberg trace formula [22].) 

Thus, the goal of the present paper is to help clarify 
the above controversy by a mainly numerical investiga- 
tion (Monte-Carlo simulation) of a non-linear cr-model in 
rf = 3 dimensions taking the HSP as coset space G/K. 
To anticipate the main result, we find a continuous phase 
transition, one-parameter scaling (IPS) and multifractal- 
ity. In other words, the data yield no evidence for a jump 
of the IPR. The localization length exponent of our model 
is obtained as i/ = 1.15 ± 0.15. 

This work is organized as follows: After a short intro- 
duction of the HSP in section II we state in section III the 
model which is to be investigated. In section IV we de- 
scribe the main ideas which enable us eventually to treat 
the problem using a Monte-Carlo (MC) algorithm. The 
results and fundamental limitations of this algorithm are 
outlined as well. We conclude with some remarks con- 
cerning the applicability of the MC algorithm for related 
problems. 



11. THE HYPERBOLIC SUPERPLANE (HSP) 

A. Geometry of the HSP 

The HSP is a supersymmetric homogeneous space 
G/K of rank one and real dimension (2,2) and can be 
viewed as a real, supersymmetric extension of the (up- 
per) hyperboloid = S0{2, 1)/S0{2). It shares several 
important aspects with Efetov's supermatrix spaces [11], 
more precisely, with the model I for the case of orthogo- 
nal symmetry. These are (besides the reality constraint) 



• perfect grading, i.e. equal numbers of bosonic (com- 
muting, even) and fermionic (anticommuting, odd) 

degrees of freedom 

• non-compactness 

• a positive (super-)Riemannian curvature 

A detailed introduction of the HSP - in particular, a pre- 
cise definition as a homogeneous space G/K- was given 
in Ref. [20] (see also Ref. [21]); the complex counterpart 
of the real HSP was introduced in [22]. These definitions 
will not be repeated here. Instead, for our purposes it 
will be sufficient to give a brief idea of the HSP and to 
put forward a useful parameterization. 

An element ip of the real HSP can be viewed as a five- 
component supervector, tp"^ = {■ipo,tpi,ip2,tp3,tp4) where 
tl^o,'4>i,'4>2 are even and tpzjipi are odd. We will use an 
involution (adjoint, complex conjugation) of the second 
kind which means that for odd elements we define 



-ii and [iiij] = ^i^j 



Using this involution we can express the reality condition 
for the HSP as 

The supervector obeys a (non-linear) constraint 
rP^Orp = 1 

where ip^ = ip and the metric 6 is given by 

^ = diag(l, -I4) . (1) 

Moreover, we require for the body ^ m(V'o) of t/Jo the con- 
dition m{ipo) > 1. We now introduce polar coordinates. 



/ cosh r 
V sinhr e 



f cos4> (1 — rjT]) \ 
sin(?!> (1 — f]r]) 

V 



J 

(2) 



with r,(f> even, < r, < < 27r and odd elements T],r]. 
Let us denote the "radial" part of a supervector ip hy , 



ij;'''^ = (cosh r, sinh r, 0, 0, 0) 
and the origin of the HSP by o, 

0^ = (1,0,0,0,0) . 



(3) 



(4) 



(The above coordinates are motivated on a group-theo- 
retical level in Ref. [20].) 



• its structure as coset space G/K 



^The body m(x) of an even element x denotes the ordinary 
part of s, which results after subtraction of all nilpotent terms 
of X. 
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in. THE MODEL 
A. Definition and motivation of the model 

In this section we introduce and motivate the model 
which is to be investigated; we also define some correla- 
tion functions of interest and specify the order parameter 
function. Finally we show how physical quantities can be 
extracted from the model. 

The model Efetov used [14] is a system of impure metal 
granules (with a large number of states) which are in con- 
tact with each other. The (macroscopic) conductivity is 
then essentially governed by the probability of the elec- 
tron jumping from one granule to a neighboring one. This 
"granular" model is closely related to Wegner's TV-orbital 
model [4]. There is some evidence [23] that the Anderson 
model and the A''-orbital model (or the granular model) 
exhibit the same critical behavior. The above model of 
disordered metal granules can be described by a super- 
symmetric non-linear cr-model on a lattice; its free energy 
(for identical granules) is written in the form [14] 

F[Q] = -jY1 ^^^QiQj -l^Y. ^^^^Qi ■ (5) 

Here, the first sum runs over nearest-neighbor granules, 
J is related to the coupling between the granules and 
the parameter (3 is given by /3 = |(w -I- ie)T:uV , where 
V denotes the volume of a granule, v the average DOS, 
e — > 0+ and w is a frequency. The supermatrices Q and 
A have dimension 8x8. Their precise forms can be found 
in [11,14]. 

The correlation function of the DOS, ifi, and the 
density-density correlator, K2, are given by 

if 1 {x, y- E) = G+{x,x;E)G-{^;^ , (6) 
K2{x,y;E;u;) = G+{x,y;E)G-{y,x;E - u;) , (7) 

where the overbar stands for the disorder average and G+ 
and G~ are the retarded and advanced Green functions. 
These correlators can be calculated via 

K,{x,y;E) oc [ JldQi Qll{x)Qll{y) e^p -F[Q], (8) 

K2ix,y;E;Lj) oc f HdQi Qllix)Qlliy) e^p -F[Q]. (9) 
i 

Here, upper indices denote 4x4 sub-matrices of the 
Q-matrices and lower indices the corresponding entries 
within these blocks. 

Now let us consider the following statistical toy model: 
At each site a; of a given lattice attach a field ^p{x) 
where is an element of the HSP. We define the action 
(resp. the free energy) of the model as follows: 



sm = 7 E i^H=^)my) - 1) + /3 E(^o(^) - 1) 

x-^y X 

= 5o[^] + 5i[V] , (10) 
P = -i{u) + ie)vV , 6^-0+ , 

where 6 was defined in eq. (1) and the notation x ^ y 
means here and henceforth that x and y are nearest 
neighbors. The "interaction term" So[ip] contains a sum 
over all nearest neighbor pairs with an inverse coupling 
constant 7 > which is a measure for the disorder of the 
system. Of central importance are the symmetry prop- 
erties of Slip]: The term Sq is invariant under a "global" 
transformation ^(a;) i-> g ■ il}{x), 5 G G . The second 
term, 5i, breaks this symmetry, but it leaves a residual 
symmetry unbroken. Namely, Si is invariant under the 
action of the stability group K C G, i.e. under global 
transformations of the form ^(x) i-> k ■ ip{x), k e K. 
(The action of elements of the stability group K leaves 
the origin (4) invariant.) Expectation values of functions 
A{ip) are defined in the usual way, 

< A{ij) > =Z~^ J [D^]M^) exp -S[iP] , 

N 

[D^] =]lD^{xi) . 

i=l 

Here, TV is the number of lattice sites, and for the parti- 
tion function Z we have 

Z = I [DrP] exp -S[rP] = exp -5[V']|^(,)=„ = 1 

because of the K-invariance which remains unbroken. 

In order to motivate the present model (10) we note 
that a detailed analysis of the supermatrix models of the 
form (5) has revealed [13,14] that the cniicaZ properties of 
these models depend crucially on their symmetry struc- 
ture with respect to the groups G and K. The toy model 
(10) on the HSP imitates exactly this symmetry struc- 
ture and is therefore expected to yield a qualitatively 
correct description of the critical behavior of the local- 
ization transition. 

The use of the HSP which has only one (non-compact) 
radial coordinate is further motivated by the observation 
that the critical behavior of the order parameter func- 
tion (see below) is governed by just this non-compact 
variable; the other (compact) radial coordinate(s) of the 
supermatrix models are of no interest in this case. This 
allows for a description of the MIT using the simplified 
toy model. (Of course, one cannot hope to use this toy 
model successfully in situations where the compact radial 
sector of the supermatrix spaces is of equal importance 
as the non-compact one.) 
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B. Correlation functions and physical quantities 

Notwithstanding the formal analogy between the mo- 
del (10) and spin models, some usual methods which are 
used in phase transition theory are not applicable here. 
These include in particular mean-field methods, which in- 
volve the magnetization, which is non-vanishing only in 
the ferromagnetic phase and identically zero in the para- 
magnetic phase. The reason is the equality < >= o 
which follows from the unbroken K-invariance so that the 
mean DOS is clearly not critical. (This is in accordance 
with the situation in the supermatrix models where the 
expectation value of Q equals A and is therefore a con- 
stant.) 

The only non-trivial 2-point function of our model is 
K{x,y;j,f3) = Kix-y;j,f3) = (-)l''l <M^)My) > 

(11) 



with 



G 1, 



ipn even 
V'm odd 



(Note that < ipo{x)ipo{y) > is not critical. This is a 
consequence of the fact that the considered model has less 
degrees of freedom than Efetov's supermatrix models.) 

The value of the coupling constant 7 decides whether 
the system is in the phase of localized or extended states 
(if both phases exist). Let 7c be the value of 7 at which 
a possible phase transition occurs. Then we have the 
correspondence 

localized states O 7 < 7c 
extended states o 7 > 7c 

Both phases differ from each other with respect to 
the symmetry breaking behavior [13]. In the disordered 
phase (insulating phase, phase of localized states) the 
broken G-invariance is restored in the limit P ^ 0. 
The conducting phase (metallic phase, phase of extended 
states), however, exhibits a scenario of spontaneous sym- 
metry breaking: the G-invariance remains broken even 
after the symmetry breaking term Si vanishes in the limit 
/3 — > 0. As was shown in Ref. [13] in the context of the 
supermatrix models this different symmetry behavior to- 
gether with the non-compactness of G/K implies a dif- 
ferent behavior of the correlation functions. One finds 



lim PK(x 

8^0 



const , localized states 
, extended states 



(12) 



The two phases of the system can be distinguished by an 
order parameter, more precisely by an order parameter 
function. [13] Let F{ip{x)) be that function which is ob- 
tained if one integrates the statistical weight exp —S[ip] 



over all sites of an (infinitely extended) lattice except the 
site X, 

F{^{x)) = JU D^iy) exp -S[^] . (13) 



For symmetry reasons F depends only on the radial part 
tp'^ (3) of tp and can be used to discriminate the two 
phases. 



lim FW) oc I / 

0^0 { B. foi 



(the constant function) , 7 < 7c 
for r — > 00 decaying function , 7 > 7c 



A physical meaning of the above function F was given 
in [18,23] where it was shown that F is closely related to 
the distribution function of the local densities of states 
(LDOS). 

For future reference let us define 



L(^,^') = exp(-7(^W-1)) , 
D0{rP) = exp(-/3(V'o - 1)) ■ 



(14) 
(15) 



These functions show the following important symme- 
tries, 

L{^,^')=L{g^,g^') V5 G G , 
D0{^P)=D0{k^P) VfcGK, 

which means that is a K-radial function, 

D^i^) = D0(r) = D0(r) . 

In addition to the correlator (11), we will also fre- 
quently be interested in "diagonal" correlators of the type 



K^J,p)=^^l,^-'Ul\x)iJT{x) 



(16) 

= {I + m)- \ l + m-2)\, 

where /, m are integers larger than or equal to zero, 
q = l + m > 2, and x is an arbitrary lattice site. The nor- 
malization factor is chosen so that the dependence of K" 
on I and m enters only via the sum q = l + m and that for 
7^-0 if approaches unity. [For q = 2 and / = 2, m = 
(or vice versa) this correlator reduces - since ipi and ^2 
are real - to the correlation function K(0; 7, /3), eq. (11).] 

Let us now demonstrate how we can (at least formally) 
extract physical quantities from the above defined corre- 
lation functions. Motivated by the unusual properties 
of the MIT on the BL we will focus our interest mainly 
onto the following quantities: the IPR (and its higher 
moments) and the diffusion coefficient D. 

Let us start from a tight-binding model with lattice 
sites x,y,. . ., normalized eigenstates ipk and an average 
level density p{E) per site at energy E. The inverse par- 
ticipation ratios (IPR's) Pq{E),q G R, 
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j d'^r\%l){r 



\2q 



are a very sensitive measure for the degree of localiza- 
tion in the system. [24] To which extent the states are 
localized, can be inferred from the dependence of Pg on 
the length L of the system (we set N = L'^, where d is 
the spatial dimension). For extended states (deep in the 
metallic region), Pq scales like 



Pn OC 7V-(«-l) (X L-(9-l)'' 



(17) 



whereas for strong disorder (which means that the lo- 
calization length ^ is much smaller than L) Pq is rather 
independent of L. In the limit of extremely localized 
states all Pg tend towards unity. Thus, the IPR's can be 
used to discriminate localized from extended states. 

Around the critical point 7c (i.e. in the range where 
^ exceeds by far any microscopic length scale Imic) one 
makes the following ansatz relying on the assumption of 
one-parameter scaling (IPS),^ 



P,(7;i)<xL--(«'/,(|), T{q) = {q-l)d*iq) (18) 



with some scaling functions fg . The "generalized dimen- 
sion (s)" d*{q) differs from the spatial dimension d and is 
a function of q. For L — > oo Pg must become independent 
of L, therefore the scaling functions fq must (for L ^ ^) 
satisfy fg{L/0 oc (L/0'"^''\ Thus, Pq(j ~ 7^; L ^ oo) oc 
^-•^(9). Using the relation ^ cx: (7^ — 7) (for 7 < 7c) 
one concludes therefore that Pg scales like 

P,(7 <jc;L^^)cK (7c - 7)"*"', Aq) = M(D, (19) 

and the multifractality of the system is reflected in a 
non-trivial dependence of 7r(g) on q. One can show in 
a rather general context (see for example [25] and ref- 
erences herein) that T{q) is a monotonically increasing 
function with negative curvature. 

The IPR's can be rewritten in terms of Green func- 
tions, [24,26] 



P2{E) = lim 



G+{x,x;E)G-{x,x;E) 



P,{E) 



v^o ■Kp{E) 

(2^)"-' (/-l)!(m-l)! 



lim 

r;->-0 



2irp 



{l + m-2)\ 



:iG+ix,x;E)y {G-{x,x;E)y 



where q = I + m and /, m G W. 

Using in a second step the techniques of supersym- 
metry, the above averaged products of Green functions 



^Of course, this IPS ansatz has to be checked explicitly by 
a numerical investigation. 



can be represented in terms of certain correlators, i.e. ex- 
pectation values of the bosonic blocks of the Q-matrices 
with respect to the generating functional F{Q), see for 
instance equations (6)- (9) and Ref. [26]. Because our 
model (10) is only a toy model there exists, however, no 
such immediate physical interpretation of the correlators 
K (11) and K''{'y;P) (16). Nevertheless, K corresponds 
formally to the correlator (9), which in turn is related 
to the density-density correlation function K2{x, y; E; u>), 
eq. (7). Further, the identity [13] 

lim T)Ki {x, y; E) = lim rjif 2 {x, y; E; 0) 

relates the correlation functions (6) and (7) in the local- 
ized regime. Therefore K plays also the role of Ki, and 
we obtain 



P2{j)= lim l3K{0;j,l3) . 



(20) 



Similarly, concerning the whole set of IPR's {q > 2), we 
can formally identify 



P,(7) = lim/39-iif''(7,/3) 

/3->0 



(21) 



which relates (at least, formally) the IPR's Pg to expec- 
tation values of the model (10). Let us note, however, 
that due to the reduced number of variables in our toy 
model (compared with the supermatrix spaces) the cor- 
relator is non-trivial only for even I and m. Thus, we 
can only hope to calculate IPR's with an even q within 
our model. 

The diffusion coefficient £'(7) can be extracted by com- 
parison with the long-range behavior of the correlator (7) 
- or, equivalently, (11) - on the conducting side of the 
transition. More precisely, its Fourier transform reads 



K{k, iv; 7) oc 



1 



D(7)fc2 



fc ^ 0, w ^ 



(22) 



from which .0(7) can be inferred. Consequently, we can 
calculate Pg(7) and .0(7) by consideration of the single 
function K. 



C. A simple illustration 

The simplest case, which can be considered, is a lattice 
with just N = 2 lattice sites, x and y, with attached fields 
V' and V'', resp. Let us evaluate the correlator (21) and 
the IPR's (21) for this particular case. We have 

Fir) = D0ir) I L{r,^')D0m ■ 

We can use polar coordinates (2) and easily integrate out 
the variables {(j),r],fj) at the site x which yields (for l,m 
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even and greater than or equal to zero) a factor uf^''] 



(;+m-2)! ■ 

In order to integrate out the remaining variables note 
that because of the infinitesimal only the asymp- 

totic domains r,r' — > oo contribute (cosh r, cosh r' ~ 
/3~^), hence we substitute 

w =/3coshr,w' = /3coshr',y>' = 4>'/l3,x' = n' 113 . (23) 

Pg(7) becomes then, using standard integral identities 
[27], 



/32 



(1 - -ywLj'x'x') 



(g-2)!^"Pr-"-V2("+2) 



(24) 



with 

Lo(w,w',v'') = exp-7 



(25) 



Note that all powers of /3 have cancelled. After a straight- 
forward calculation one obtains for the last remaining 
integral 

where (j) is the error function. In particular, for q = 2 the 
result is 

^2(7) = 1-^^62^ (1-<^(V^)) ■ 

The limiting behavior of Pq{'y) for 7 — > and 7 — > 00 
can easily be read off. 



{ 2-(5-i) \ 



7^0 
7 — > 00 



as required in eq. (17). 

Let us now turn our attention to the three-dimensional 
case where integrals of the above type can, of course, be 
solved only numerically. 



model (10) numerically. For simplicity, we will restrict 
ourselves to diagonal correlation functions; specifically 
we will calculate the set of IPR's -Pg(7) (21). In princi- 
ple the algorithm is capable of calculating non-diagonal 
correlation functions as well. 

In order to apply MC techniques it is absolutely nec- 
essary, however, to find a method which treats the fermi- 
onic variables appropriately. It is well-known that at- 
tempts to study systems which include dynamical fermi- 
ons by means of MC simulations face severe problems. 
This is mainly due to the so-called "minus sign" problem 
which arises from the presence of Grassmann variables 
that make a definition of a positive-definite probability 
measure problematic [28]. However, it turns out that an 
MC simulation of the model (10) is possible in spite of 
the presence of fermionic variables, because of the self- 
terminating property of the Grassmann polynomials. 

The underlying idea of the algorithm has been sug- 
gested some years ago [13] but has never been pursued. 
Thus, no numerical simulations of the kind below have 
yet been performed. 



A. Derivation of the Algorithm 

To begin with, let x be an arbitrary lattice site with 
attached field ^. We write the IPR P,(7) (21) in the 
following form 

P,(7) =M,m-Mim/3''-iy D^^'^^Fir) , 

Fir) = D0{r) I n my)D0{r{y)) (26) 

where q = I + m. (Pqij) does not depend on the specific 
choice of x, of course.) Similar as in section IIIC we 
integrate out the angular and fermionic variables at the 
site X. The substitution w = /3coshr leads afterwards to 



PoH) = lim / dLj(x)( 

0^0 J 



U)(x) 



9-2 



(a -2)! 



n L{r{x),rP{x')) n L{rP{xi),rP{xj)) 




yjix 



IV. MONTE-CARLO SIMULATION 

^In the limit /3 — )■ no boundary terms arise when one 
In this section we introduce a MC algorithm which changes to polar coordinates. - The ^-integrations run from 
allows us to calculate some correlation functions of the to cx3. 
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with 



{{...)) = I Mx)e--^'^i{r{x)) (...), 

I{V(x)) = I YlD^Piy) Y[ L{r{x)Mx')) 

Note that the normalization of this average is 

Z = {{!)) = j da;(a;)e-'^(^)/(^''(a;)) = 1 

because I{ij}^{x)) = 1 due to the absence of symmetry 
breaking terms. 

Thus, the correlation functions can be expressed as 
"thermodynamic" expectation values, similar to the con- 
cept of statistical mechanics. Thus, the following scheme 
becomes apparent: 

• Generate configurations according to a certain dis- 
tribution function which will be specified below. 

• Average the symmetry breaking terms 



{q-2)\ 



over these configurations. 

Next comes the problem of how to deal with the Grass- 
mann variables. Performing at each lattice site the sub- 
stitution (compare (23)) 

{r,(f),fj,r)) ^ (w = /3coshr,y> = (j)/0,x = 'n/0,x = v/0) 

we integrate out all remaining Grassmann variables 
{Xi x) '■ First note that as a direct generalization of 
eq. (25) we find 



and ibj < ^'^^ occupied 1 "interaction term" ac- 

■' [ occupied J 



(27) 

where Axij = Xi~Xj- Now let us call a bond between ipi 
on-occupii 
occupied 

cording to eq. (27) is given by < . 

Then, after integrating out the Grassmann variables, one 
can express I(ip'^{x)) as a sum over all possible graphs 
(consisting of occupied and unoccupied bonds) on the 
underlying lattice. The following simple graphical rules 
for the evaluation of I{il)^{x)) hold: 

1. Each graph which contains a loop of occupied bonds, 
vanishes. The reason for this is the identity 

AX12AX23 • --AXnl = . 

Thus, only tree graphs contribute to J(t/;''(a;)) . 



2. Each graph that is not linked to the site x (via a 
cluster of occupied bonds ) vanishes. The reason for 
this is that for a non-vanishing contribution the 
number of Grassmann variables (provided by the 
occupied bonds) must match exactly the number 
of fermionic integrations (provided by the vertices 
except by the site x). 

Thus, only tree graphs that emerge out of the site x, can 
contribute. Let us denote this set of contributing tree 
configurations by {G}. 

After we have integrated out the Grassmann variables, 
the partition sum can be rewritten as 



Z = Y.f d^{x)e-^(<^^ [ n 

{G}-' •' yi^x 



{G}' 



27r 



e{G} 

JJ Lo{wi,LLjj,(pi - (Pj) Y[ l^ki^l- 

Xi^^^Xj X]^^^^X\ 

In order to circumvent the - from a numerical point of 
view - inconvenient condition w > together with the 
singularity of the measure at a; = the following substi- 
tution appears useful, 

\{x) = exp— w(a;) , t{y) = lnw(2/) (for y ^ x). 

Thus, we have 

A(7) = E / mM4>)v{(i>,G), 

df{y) 

V{4>,G) = 7^-iw(a;)^(") JJ a;(2/)''^^"' 

y^x 

where z{xk) is the local coordination number of the site 
k, i.e. the number of occupied bonds at k. 

In principle one could now calculate Pqij) via an MC 
algorithm on the configuration space {(wj, (p,), {G}} . 
However, such simulations have revealed that the relax- 
ation into equilibrium is - even upon small lattices - very 
slow, and far from being satisfactory. It is therefore ex- 
pedient to reduce the number of degrees of freedom of 
the model further, by integrating out the angular vari- 
ables ifi. This integration is just an {N — l)-dimensional 
Gaussian integral with the additional constraint (p{x) = 
and yields 



/»00 PC 

D[cj>] = / rf(e-'^(^)) n / d{\nu;{y)) / 

Jo „,-L„ J —00 J —{ 



n 

k^x 



d(pk 
2tt 
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Here, the matrix results from the real, symmetric, 
positive definite matrix A with 







z = J 
i ^ j 
otherwise 



(28) 



by erasing the line and column corresponding to 
the lattice site x. (Note that det^j; > 0.) 
Therefore we can reduce the configuration space to 
{[Q] = {wi, . . . , wjv}, {G}), and obtain finally 



n(7) = E / D[m,mwin,G) 

{G} 

D[n]= / d(e-'^(^))n / d{\na;{y)) 

Jo „^^J-oo 



A,iQ) 



--— llexp-^(2/) 



N-l 

W{n,G) = (^) ' {detA^y-^ivix) 



z{x) 



L^{ui,u}j) = exp 



.27r 

2 U!i UJi 



(29) 
(30) 
(31) 

(32) 
(33) 



Notice that Y.{g} I D[^]W{n,G) = 1 and that the ap- 
pearance of the exponents z(y) — 2 is "natural" because 
the average (local) coordination number z of a tree con- 
figuration is given by 

The scheme of a MC simulation is now (see for ex- 
ample Ref. [29]) to approximate the integral (29) by its 
importance sampled average 



_ 1 ^ 



(34) 



where the M configurations {Qi,Gi) are distributed ac- 
cording to their statistical weight, i.e. with probability 
Pini,G,)o:W{n„G,) . 

The above considerations motivate the following MC 
algorithm: 

1. Choose an arbitrary lattice, a site x and a start 
configuration ([fi], G) of fields and bonds (the latter 
exhibiting a tree structure). 



2. Generate (locally, i.e. by changing the value of a 
single, randomly chosen w) a new field configura- 
tion ([n'],G). Calculate the weight ratio ^ 
and accept the new configuration with probability 
Psite =mm(l, ty^n.G) )■ 

3. Choose randomly a nearest-neighbor pair (ri,r2) 
where ri and r2 are not linked with each other by 
an occupied bond. Generate a new tree configura- 
tion G' by breaking the bond between r2 and its 
predecessor and linking ri and r2 if r2 is not an 
"ancestor" of ri , which means that the path from 
ri towards x along occupied bonds does not visit 
r2. This procedure is illustrated in figure 1; it guar- 
antees that the new bond configuration is again of 
tree type. Calculate the weight ratio of the new 
and the old bond configuration, 

w(ri) 



w{n,G') 

VF(fi, G) a; (predecessor of 



and accept the new configuration with probabil- 
ity pbond = min(l, ^^(q' )• Note that all possible 
tree configurations can be attained within this pro- 
cedure. 



r2 



rl 












> 











r2< 
I 
I 
\ 



r2^ 
/ 

I 

\ 



-^ 








> 




— 








X 








> 






; 
I 







FIG. 1. Example of contributing bond configurations 
which represent trees emerging from the site x. Left: r2 is 
an ancestor of ri; right: r2 is no ancestor of r\, therefore a 
new tree configuration evolves. 

Before I present the results of numerical simulations, it 
seems reasonable to get first a qualitative understanding 
of the model (29)- (33). I will discuss this point taking 
the one-dimensional chain as an example. 
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B. Qualitative Discussion for a Id chain 

In a one-dimensional chain with N lattice sites the 
sum over trees reduces to a triviality since only one tree, 
namely the chain itself, contributes. We will therefore 
omit the index {G} in the sequel. 

Consider first the limit 7 — > oo . In this case the factors 
L^{uJi, ujj) in the weight W{Q) force all uJi to take roughly 
the same value which I denote by u>. Consequently one 
can approximate 

and up to factors which are independent of u) one finds 
W{n) const X w^-i JJ 6{uJi - Uj). 

In the evaluation of the quantity Pg all integrals ex- 
cept the integral over w(a;) break down due to the delta- 
function and one arrives at 

dQexp{-Nu^)j—-^^ = 7V-(«-i), 

which is the correct limiting behavior. 

For finite values of 7 the statistical mechanics of the 
model is governed by two competing effects: On the one 
hand, the variable liJx takes values of order unity, since 
/-D[n]W(f2) contains the factor a;(a;)^ exp — a;(a;). On 
the other hand, the radial coordinates at the boundaries 
of the chain, uii and wjv, say, show a strong singularity 
of type 1/w in eq. (32) and tend therefore to zero. The 
inverse square root of det does not resolve this singu- 
larity. Depending on the value of 7 the intermediate sites 
"feel" the effect of the dynamics of (v^, Ui and wjv and 
balance it. For small 7 the effect of the "fixing" of the 
site X is only weak, and consequently all - except u)x 
- are pushed to zero. Thus, all IPR's Pq tend to unity. 

For lattices which are of non-tree type (for instance hy- 
percubic lattices m d>2) the qualitative picture is simi- 
lar; in addition to the previous discussion any site i with 
an (at the moment) large value of uji (compared to the 
other lattice sites) will accumulate occupied bonds at the 
cost of those sites j with a small value of ujj. Conversely, 
a large number of occupied bonds at a given lattice site 
i will favor a large value of . 

C. The Calculation of the Determinant 



the one-dimensional chain with N sites it can be calcu- 
lated with 0{N) multiplications. However, for a hyper- 
cubic lattice of size N = L x L x L and skew-symmetric 
(i.e. helical) boundary conditions ^, a naive implementa- 
tion using standard libraries would take 0{N^) = O(L^) 
steps, which is clearly unacceptable. The reason is that 
although the matrix A itself is sparse, it becomes in gen- 
eral dense after an LCZ-decomposition (which is used for 
the calculation of its determinant). Such a decomposi- 
tion generates O(iV^) "fills" in the matrix, i.e. entries 
that were originally zero become non-zero after the LU- 
decomposition. 

In a first improvement, one can modify the boundary 
conditions in such a way that the matrix A displays a 
band structure with bandwidth 0{L^). This modifica- 
tion affects only a fraction of 0{N^3) lattice sites - for 
which the boundary conditions are changed from periodic 
to free (hard) - and is therefore negligible in the thermo- 
dynamic limit. We can then use a standard Cholesky 
factorization for the determinant of a real, symmetric, 
positive definite, banded matrix. But since the CPU-time 
in this case is still of 0{L^) such algorithms are practi- 
cable only for very small sj^stems (L < 5). 

One can, however, calculate the determinant of A^ in 
an approximate but controlled manner in 0{L^) time. 
First, one splits the matrix A (28) into its diagonal and 
off-diagonal parts, 

A = B-C,Bij = ar^Sij,Cij = Si^j = { J ' otherwise " 

For each matrix T, let denote the matrix which results 
from erasing the line and column x. Using the identity 
det = exptrln, and expanding the logarithm, one then 
has 

det A^ = det det(l - B~^C^) 

= det Bx exp ^- ^ ^t^M^^ , (35) 

Mij = y/ajaj Si^j 

It remains to calculate trM^ for n = 1,2,... Note that 
M3; does not display a simple structure, due to the effect 
of erasing one line and column. However, if we redefine 

Mij = ^Jaittj Si^j , a^; := (36) 

(i.e. the line/column which is to be erased is replaced by a 
zero line/column) then it holds true that trM" = trM". 




This section deals with the "bottleneck" of the algo- 
rithm, namely the (technical, but numerically extremely 
important) question how to evaluate efficiently the deter- 
minant det^a;, which appears in W{il,G), eq. (32). For 



''By skew-symmetric boundary conditions we mean that the 
lattice sites are numbered from to JV — 1 where for each site 
k its nearest neighbors are located at the positions k±l, k±L 
and k±L^ (modulo L^). 
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Thus, we need not bother about the somewhat unpleas- 
ant structure of M^; instead, the problem of calculating 
det Ax is essentially reduced to the evaluation of 

trM" = ^M;j= J2 aiak^ak2...ak„_^ 

for positive n. The latter sum can be interpreted as a sum 
over all closed paths of length n starting at (and returning 
to) a site i on the underlying hypercubic lattice. Having 
this picture in mind, one clearly recognizes that 

trM^""^ = , 71 = 1,2,.... 

In addition, it can be shown that the series (35) converges 
at least like -2~", which allows one to truncate the series 
after a finite value of n. In practice, 2n = 4 is sufficient 
for an accuracy of better than 1%. Consequently, we can 
express the determinant of by 

det ~ [ Yl I exp - QtrM^ + ^trM^^ 

Mij = ^aittj 6i^j, tti = , = 0. (37) 

Moreover, the remaining error largely cancels itself out 
when one calculates the ratio ^^qQ§^, so that we end 
up with a relative accuracy of about 10^^ . (In order to 
achieve a performance of 0{L^) essential multiplications 
in the calculation of det j4;i;, it is, of course, necessary to 
use the sparsity of the matrix M in the implementation 
of the matrix multiplications in eq. (37). 

D. Results of the Monte-Carlo Simulation 

We have now gained a qualitative understanding of 
how the model works and how a numerical implemen- 
tation is in principle done. I will therefore present the 
results of the above algorithm. 

The algorithm has been checked at the example of a 
one-dimensional chain with N sites, since exact results 
for Pqij) are available for both N = 2 and N = oo. 
The agreement between analytical and numerical results 
is very good. 

From the analytical calculations of the model (10) on 
the BL [19] we expect the critical disorder jc for the 
3d hypercubic lattice in the range between 10~^ and 
10~^. We focus our attention therefore on this order of 
magnitude of 7 and calculate Pq(-y;L) for system sizes 
L = 4,5,6,7,8,10. ^ This is done as follows. To 



obtain a set of Pq{-y;L) (for fixed 7 and L, but vari- 
able = 2, 4, 6, . . .) we perform typically some 10^ — 10'' 
Monte-Carlo steps (MCS), where each step represents a 
complete sweep i.e. consists oi N = random updates 
of sites and bonds as described in section IV A. We 
store the pair {ijj{x),ex'p — '^y^^uj{y)) after each tenth 
MCS. From this data one can afterwards average (34) 
the expression Aq{fl) (31) for each q separately and ob- 
tains thus Pg(7; L) for all q. In order to make sure that 
the above data-points Aq{Q.i) (where i labels the MCS), 
which are to be averaged, are statistically independent, 
we calculate the autocorrelation function 

C{ti) = Aq{ni)Ag{ni+„) - Tq (38) 

and take for averaging only each {NcoTrf^ data-point 
Aq(^i). Here, Ncorr is defined as the length on which 
C{n) has collapsed to a relative value of 10%. Given 
these Neff = Ntotail^corr data-poiuts, the statistical er- 

ror scales like ^^.f^f 

After we have produced for each q and L the "raw 
data" yi = Pq{'ji;L) (for discrete 7^) together with their 
statistical errors o-j, this data is subjected to a x^-fit 
procedure. This procedure yields a continuous function 
Pq (7; L) that interpolates between the values by a least 
squares fit (see for instance [30] or an appropriate text- 
book). As an illustration figure 2 shows the resulting set 
of functions for q = 2 and L = 4, 5, 6, 7, 8, 10. 




0.01 0.1 

gamma 

FIG. 2. The "raw data" P2(7i;i) together with their sta- 
tistical errors, and the interpolating function P2{'y;L) as cal- 
culated from a x^-^*- System sizes are L = 4,5,6,7,8,10 
(from top to bottom). 

These functions tend in fact to the correct limits for 
7 — > and 7 00 and do not intersect each other. 



E. Multifractal Analysis 

^It turns out to be somewhat disadvantageous that the pre- 
cise value of 7. is not yet known which hampers a direct in- The goal of this section is to investigate the expected 
vestigation of the critical regime right from the beginning. multifractal structure of the MIT, by extracting from the 



10 



set of functions Pg (7; L) the generalized dimensions d* (q) 
(18) and the correlation length exponent u via a finite size 
analysis. 

We determine jc and the set T{q) (from which d*(q) fol- 
lows immediately) by multiplying each function Pq{'y; L) 
with a factor - with a "guess" f (g). If (and only if) 
f(g) = r(g) then all curves L'"(5)Pg(7; L) (for fixed q and 
variable L) intersect in one single point at a nontrivial 
value 7 = 7c. This is due to the expected IPS behavior 
in the critical regime (where ^ ^ Imic), 



(39) 



compare eq. (18). (T(g) does not depend on 7.) Usu- 
ally, the above "fitting" procedure is not very accurate, 
since the critical coupling constant 7c is not known and 
therefore two parameters {T{q) and 7c) are determined 
simultaneously, which might yield quite a large uncer- 
tainty. Note, however, that 7c must take the same value 
for all q, therefore this approach allows for a rather pre- 
cise estimation of 7c. Thus, for each q, effectively only 
one parameter, T(q), is fitted. (Once the value of 7c is 
known, one could alternatively extract the set T{q) from 
a plot of lnPg(7c; L) versus InL, 

In Pg (7c; L) = —T(q) In L -I- const 

which is a straight line with (negative) slope T{q).) 

Proceeding along this scheme we observe indeed IPS 
and multifractality at the critical point: the dimensions 
d*{q) = t{q)/{q — 1) are not constant but fall off with 
increasing q, approaching a saturation value doo — 0.3. 
In figure 3 we show as an example the IPS behavior for 
q = 4. Here, we have multiplied the function P4(7, L) 
with a factor t(4) ~ 1.8 which results (roughly) in a 
single intersection point around 7c — 0.04. 



0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 

gamm 



FIG. 3. The curves L"(«'Pg(7;L) for q 
L = 4, 5, 6, 7, 8, 10 and r(4) = 1.8. 



where for different q the numerically obtained values for 
7c fluctuate a little bit around that value. The set T{q) 

is displayed in flgure 4. 




FIG. 4. The obtained values for r(g) 
with q = 2,4,6,8,10,12,14. The trivial results r(G) = -d 
and t(1) = are added. 

In principle we can from this data determine the so- 
called /(a) spectrum, [25] 



/(a) = a{q)q - T{q) , a{q) 



dTjq) 
dq 



(40) 



In our approach, however, q is not a continuous vari- 
able, as we do not calculate the quantum mechanical 
wave function (from which Pg for any value of q could 
be deduced). Instead, q takes only discrete (and even) 
values. We therefore omit the calculation of the /(a) 
spectrum, since the regime a ~ qq (the maximum of 
/(a)) cannot be probed by our approach. 

The localization length exponent u can be obtained 
if one makes a guess D and plots Pg(7; L)L'''') versus 
L^/" for 7 ~ 7c and fixed q. For, from eq. (39) one 
recognizes that, if i> = i/, then the curves Pq{'j; L)L^^''^ 
have for different L (but fixed q) the same dependence on 
7— 7c and fall together, as long as the criticality condition 
^ ^ Imic is satisfied. 

We obtain 



1/ = 1.15 ±0.15 



compare figure 5. 



We find 



7c = 0.038 ± 0.004 . 
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0.1 I ^ 1 ^ ^ 1 

-0.1 -0.05 0.05 

(gamrna-gamma_c) times L'*{1/nu} 

FIG. 5. L"'''P2(7;i) vs. (7 - 7e)Li/'' ^^^2) = 1.45, 
v = 1.15) for different L. 

Although this value of (and the findings for d*(q)) 
is surely not the result of a high-precision calculation, it 
fits well into the "landscape" of estimates for v for the 
3d Anderson model that have been obtained in the past 
by means of several different methods, see for instance 
[31-35]. These values range - roughly - from v ~ 0.97 
to u ~ 1.5 together with some quite strongly deviating 
estimates. It is, of course, not guaranteed that the local- 
ization exponents of our toy model and of the Anderson 
model coincide exactly. But the qualitative agreement of 
its numerical values is a further indication that the above 
toy model yields a qualitatively correct description of the 
localization transition. 



More problematic seems the structure of the model 
at the opposite limit, that is, for 7 > 7c and large L. 
As the symmetry breaking term Ag (31), which is to be 
averaged, is essentially of the form exp — uji, the fluc- 
tuations of this quantity become enormously large if L 
is large and the are not small, which means that 7 is 
not small. Therefore one needs a large number of MCS, 
since - roughly speaking - the distribution of Ag is es- 
sentially of log-normal type, and thus the less probable 
values dominate its average. 

One might ask whether this problem cannot be circum- 
vented if - speaking in general words - instead of aver- 
aging a strongly fluctuating quantity Ai according to a 
weight distribution Wi one averages, say, ^/Ai according 

to WiyfAi, 

and hopes that the fluctuations are damped by the square 
root. However, numerical tests have shown, that this 
procedure yields no improvement. Instead, the problem 
of exponentially large fluctuations seems to be inherent 
to the MIT (at least in the formulation of a non-linear 
(T-model). 

It is the last restriction that hampers me from inves- 
tigating systems with L > 13, say. However, as we have 
seen, the IPS sets in already for L < 10 so that the inves- 
tigation of still larger systems seems to be unnecessary. 



F. Limitations of the Algorithm 

Let us briefly discuss the limitations of the applicabil- 
ity of the above outlined algorithm. A natural limitation 
is given by a maximal system size L < 15, because of re- 
strictions in CPU time. Further, because of the "critical 
slowing down" one needs more and more MCS to investi- 
gate the critical point. Apart from that, there are other 
restrictions which are specific to our model. 

First, for 7 — )• 0, all tend to zero, but non-uniformly. 
Thus, the exponent of Ly (33) fluctuates strongly - 
in spite of the smallness of 7 - due to the "almost- 
singularity" which appears in Wj/wj. Consequently, the 
weight ratio W(Q', G) /W{fl, G) might take very small or 
large values and produce a numerical underflow or over- 
flow. This problem can be resolved if one chooses in each 
MCS only small deviations from the previous conflgura- 
tion so that the large numbers cancel each other mostly. 
This is, of course, at the cost of CPU time, since one 
needs more MCS to "scan" the configuration space. For- 
tunately, these "precursors" of the singularity at 7 = 
do not affect the present analysis too seriously, since 7c 
is still large enough so that simulations around 7 ~ 7c 
are not plagued by this feature. 



V. DISCUSSION 

In this work, we have put forward a non-linear cr-toy- 
model which we assume to give a qualitatively correct 
description of the localization transition. After a short 
discussion of a simple case, in which this model is ex- 
actly solvable, we outlined an MC-algorithm and pro- 
vided numerical data which demonstrate that the local- 
ization transition within the supersymmetric formulation 
is a continuous phase transition which can be described 
by IPS. The multifractal nature of the transition has 
been confirmed. 

It would be an interesting future task to transfer the 
MC algorithm from the toy model of the present work 
onto the supermatrix models, for instance onto the uni- 
tary model. It seems that the additional degrees of free- 
dom, that arise in such a transfer, are not too difiicult 
to handle.^ The essential features of these models are 



However, due to the presence of more than one Grass- 
mann variables together with their conjugates, the contribut- 
ing bond configurations in section IV are no more necessarily 
of tree type. 
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already encountered in the present MC algorithm. Such 
a transfer would be advantageous in the sense that a 
direct connection between physical quantities and corre- 
lation functions of the model is available. Moreover, one 
can calculate the IPR's Pg for all (positive integer) q not 
only for q even. 

One could also utilize the above MC algorithm in a 
very similar way in order to calculate the order parame- 
ter function F(ip''), eq. (26) in d = 3. This yields directly 
the distribution function of the local amplitudes of wave 
functions [36,37]) instead of each of its moments. Fur- 
ther, one can investigate the case d = 2 (see for example 
Ref. [38]) where it is possible to simulate systems of larger 
(linear) size than in three dimensions. 
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the German National Scholarship Foundation for finan- 
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